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ABSTRACT 

Active galactic nuclei (AGNs) and quasars are important astrophysical objects to understand. Re¬ 
cently, microlensing observations have constrained the size of the quasar X-ray emission region to be 
of the order of 10 gravitational radii of the central supermassive black hole. For distances within a 
few gravitational radii, light paths are strongly bent by the strong gravity field of the central black 
hole. If the central black hole has nonzero angular momentum (spin), a photon’s polarization plane 
will be rotated by the gravitational Faraday effect. The observed X-ray flux and polarization will 
then be influenced significantly by the strong gravity field near the source. Consequently, linear grav¬ 
itational lensing theory is inadequate for such extreme circumstances. We present simple algorithms 
computing strong lensing effects of Kerr black holes, including effects on polarization. Our algorithms 
are realized in a program “KERTAP” in two versions: MATLAB and Python. The key ingredients 
of KERTAP are: a graphic user interface, a backward ray-tracing algorithm, a polarization propaga¬ 
tor dealing with gravitational Faraday rotation, and algorithms computing observables such as flux 
magnihcation and polarization angles. Our algorithms can be easily realized in other programming 
languages such as FORTRAN, C, and C-I--I-. The MATLAB version of KERTAP is parallelized using 
the MATLAB Parallel Computing Toolbox and the Distributed Computing Server. The Python code 
was sped up using Cython and supports full implementation of MPI using ‘mpidpy’ package. As an 
example, we investigate the inclination angle dependence of the observed polarization and the strong 
lensing magnification of AGN X-ray emission. We conclude that it is possible to perform complex 
numerical-relativity-related computations using interpreted languages such as MATLAB and Python. 
Subject headings: accretion, accretion disks — polarization — gravitational lensing: strong — quasars: 
supermassive black holes 


1. INTRODUCTION 

Accretion disks are one of the most spectacular phe¬ 
nomena in modern astrophysics. The large luminosity 
of active galactic nuclei (AGN) and quasars (as large as 
10^^ erg/s) is believed to be the result of gas accreted 
by central supermassive black holes. Up to ~I0% of the 
accreting mass can be emitted as radiation during the 
accretion process (much higher than nuclear fusion pro¬ 
cesses, about ~0.5%). Besides emission in the infrared, 
optical, and UV bands, most AGN emit X-rays. Un¬ 
like optical emission, which can be explained by stan¬ 
dard accretion disk theory (Shakura & Sunyaev 1973), 
the physical mechanism of AGN X-ray emission remains 
an enigma, despite decades of research effort (e.g., Krolik 
2008). According to the standard theory, the tempera¬ 
ture of an AGN accretion disk is not hot enough to emit 
X-rays. The X-ray emission is thought to be generated 
through inverse Gompton scattering of disk photons with 
hot electrons in the so-called X-ray “corona.” However, 
both the physics and geometry of this mysterious X-ray 
corona are not well understood. Observations have re¬ 
vealed that photons of different frequencies correspond 
to emission regions of different sizes: the higher the fre¬ 
quency, the smaller the emission size and the closer it 
is to the central black hole. Very recently, observations 
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in quasar microlensing (Kochanek 2004; Blackburne et 
al. 2006, 2014, 2015; Morgan et al. 2008, 2012; Ghar- 
tas et al. 2009; Dai et al. 2010; Ghen et al. 2011, 2012; 
Mosquera et al. 2013; MacLeod et al. 2015) have con¬ 
strained the AGN X-ray emission size to be of order 
10 gravitational radii {vg = GM/c^ where M is mass 
of the central black hole). Within distances of a few 
Tg, Einstein’s general relativistic gravity theory plays a 
fundamental role in understanding accreting phenomena. 
In particular, the X-ray emission is strongly lensed by 
the central black hole before it escapes the gravity field 
and arrives at a distant observer. It is well known that 
photons follow null geodesics in curved spacetimes. The 
spacetime around a rotating black hole is described by 
the Kerr metric (Kerr 1963). The strong field environ¬ 
ment of AGN X-ray emission invalidates the linear ap¬ 
proximation as made in standard gravitational lensing 
theory (Schneider et al. 1992). Consequently, numeri¬ 
cal integration of the geodesic equations is necessary for 
many problems involving the Kerr metric. Another well- 
known fact is that a photon’s polarization vector is paral¬ 
lel transported along the photon’s geodesic and its plane 
of polarization is rotated as the photon passes the rotat¬ 
ing black hole, i.e., the so-called gravitational Faraday 
rotation (Balazs 1958; Plebanski 1960; Ishihara et al. 
1988; Agol 1997; Agol & Krolik 2000; Frolov & Shoom 
2012; Yoo 2012). Since X-rays are emitted in regions very 
close to the black hole, the observed X-ray polarization 
will be influenced by the strong gravity field. Following 
the pioneering work of Stark, Connors, and Piran (Con- 
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nors & Stark 1977; Stark & Connors 1977; Connors et al. 
1980), X-ray polarization has now become a very promis¬ 
ing probe of important parameters such as the black hole 
spin and the inclination angle of the accretion disk (Agol 
1997; Agol & Krolik 2000; Dovciak et al. 2004; Li et 
al. 2009; Schnittman & Krolik 2009; Schnittman et al. 
2013). It is therefore essential to include the effects of 
strong gravity when studying AGN X-ray emission. 

Numerical ray-tracing in the Kerr metric is not new 
(Cunningham 1975; Luminet 1979; Rauch & Blandford 
1994; Bromley et al. 1997; Beckwith & Done 2004; Fuerst 
& Wu 2004; Schnittman & Bertschinger 2004; Broderick 
2006; Cadeau et al. 2007; Dolence et al. 2009; Psaltis & 
Johannsen 2012; Chan et al. 2013; Chen et al. 2013a, b; 
Schnittman & Krolik 2013). However, public codes are 
relatively rare. Exceptions are Dexter & Agol (2009; a 
FORTRAN code, fast but mathematically convolved; see 
also Yang & Wang 2013), Vincent et al. (2011; a C-|—|- 
code, not parallelized), Kuchelmeister et al. (2012) and 
Chan et al. (2013; both are GPU-based ray-tracing code 
solving 2nd order geodesic equations using GUDA). Po¬ 
larization and gravitational Faraday rotation have been 
studied extensively in the literature (Agol 1997; Agol & 
Krolik 2000; Schnittman & Krolik 2009); however, there 
seems to be no public code with a focus on X-ray polar¬ 
izations. Furthermore, there seems to be no ray-tracing 
code written in interpreted high-level languages such as 
MATLAB or Python. We present simple algorithms for 
computing the effect of strong lensing of AGN accretion 
disks by a Kerr black hole, including effects on polariza¬ 
tion. Our algorithms are realized in the high-level lan¬ 
guages MATLAB and Python, and they are fully paral¬ 
lelized and can be easily realized in other languages such 
as FORTRAN and C-b-b. 

The plan of this paper is as follows. An introduc¬ 
tion to the Kerr metric, the tetrad formalism, the strong 
lensing formalism, and the polarization formalism is 
given in Section [51 The graphic user interface (GUI) 
and numerical algorithms are discussed in detail in Sec¬ 
tion |3l In Section |4l we give an example of the use of 
our MATLAB/Python program KERTAP (KErr Ray- 
Tracing And Polarization). In Section [5] we conclude. 

2. THEORETICAL FRAMEWORK 

Since the gas accreted by black holes carries angular 
momentum, the black holes powering accretion processes 
are expected to possess angular momentum themselves. 
The proper metric describing the spacetime around a ro¬ 
tating black hole was discovered by Kerr in 1963 (Kerr 
1963) and named after him. We introduce the Kerr met¬ 
ric in this section, and discuss the other formalism needed 
in the strong lensing and polarization analysis of AGN 
X-ray emission. 


2.1. Kerr Metric 

In the Boyer-Lindquist coordinates, the Kerr metric 
can be written as 

y;2 2 

ds^ = —a^dt^ H- TT sin^ 6{d(j) — Lodt^ + -^dr^ -b p^dSft) 

p‘^ A 

where (t, r, 0, (jj) are the four independent coordinate vari¬ 
ables, and we have defined 

p =r a COS 0 


A = r^ 

S2 = p2 


-2Mr + a^ 
(r^ + a^) + 





S2 

2aMr 


E2 


2Mr 



( 2 ) 


with constants M and a the mass and angular momen¬ 
tum (i.e.,J/M) of the Kerr black hole (we have taken 
c = G = 1). The geodesic motion of a photon (its 4- 
momentum satisfies = 0, the index a runs from 1 to 
4, with 1 the time component) can be described by the 
group of 8 first order Hamilton equations (Arnold 1978) 

(ia;“ _ dn 

d^ dpa 

^ = ("U 

d^ dx<^' ^ ’ 


where ^ is the affine parameter and 'H{x°‘,pb) = 

\g°‘^{x)paPb is the Hamiltonian. The number of equa¬ 
tions to integrate reduces to six since pt and are mo¬ 
tion constants (the Hamiltonian % does not depend on 
t ov (p). There is another motion constant discovered by 
Carter (1968), 


C =pI + cos^ 9 




( 4 ) 


The number of ordinary differential equations (ODEs) 
used to solve the geodesic equations could have been 
reduced by one using this constant; however, we use 
C as an independent check of our ray-tracing code. A 
Kerr black hole has two horizons (outer and inner) de¬ 
fined as r± = M ± \J M2 — o?. The maximum angular 
mome ntum of a Kerr black hole was found in iThorne 1 
(fT^ to be Umax = 0.998. We stop the ray-tracing 
when the geodesic is within 0.1% of the outer horizon r_|_. 
For accretion disks (normally assumed in the equatorial 
plane, 9 = 7r/2), an important quantity is the so called 
innermost-stable-circular-orbit (ISCO) in the equatorial 
plane (Chandrasekhar 1983). For example, nsco = 
and 1.24 Tg for a Schwarzschild black hole (a = 0) and 
an extreme Kerr black hole. We use risco as the inner 
cutoff of the accretion disk in this paper. 


2.2. Tetrad Formalism 

In curved spacetime the frequency (or energy) of a 
photon depends on the spacetime point where the fre¬ 
quency is measured, and also on the 4-velocity of the 
observer measuring the frequency {u°'Ua = — 1). If Pa 
is the 4-momentum of the photon, the observed photon 
frequency v = —u°‘pa. Let Uo and Vf. be the photon’s fre¬ 
quency measured by a distant observer and an observer 
comoving with the light source (emitter). We define the 
redshift factor 

Vo {-U°‘Pa)o 

9= — = -f - (5) 

Ve (-lt“Po)e 

where and are the 4-velocity of the distant observer 

and the source. 

For the distant observer, we choose u“ = i(l,0,0,a;), 
the so-called zero-angular-momentum observer (ZAMO; 
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Bardeen et al. 1972) where w and a are defined in Eq. 

A photon’s momentum can be specified using the or¬ 
thonormal tetrad of ZAMO instead of the coordinate ba¬ 
sis, i.e., 


0 
0 

0 0 p 0 

S sin 6 Q Q X) sin 6 


a 0 0 

0 ^0 
VA 


\ fP"' 


P 


( 6 ) 


We will use Eq. (jH]) to compute the polarization measured 
by a distant observer (see Section 1^31) . The four-velocity 
of the source is model-dependent. Eor accretion disks, we 
assume Keplerian flow for the accreting gas (Bardeen et 
al. 1972), i.e., 

= 7(1, 0,0,11^), (7) 


where 


drop the factor p for sources at cosmological redshifts be¬ 
cause these sources are observed with tiny solid angles. 
If we ignore Kerr strong lensing and Doppler shifts, pho¬ 
tons travel along straight lines, and there is no frequency 
redshift. Consequently is a conserved quantity along 
the light path, 

^unlensed ^ f ^ / /:----(n)d4l2) 

J ’’obs J 

where n is the 3D photon direction at the observer, dA is 
the differential area element of the accretion disk, and we 
have changed the solid angle integration at the observer 
into a 2-D surface integral over the accretion disk. On 
the other hand, for Kerr lensing, we have 

F“ = J gHn)Itr^%x^ih),pain))dno, (13) 


o - 

K= r-i/2 A 

7= (-ffii - (8) 


We similarly define an orthonormal tetrad in the comov¬ 
ing frame of the Keplerian flow, and to change between 
tetrad and coordinate components we use 


/p‘ \ 


p^ 

- 


\p^ / 



where 
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(9) 


( 10 ) 


The above equations are useful since both the intensity 
and the polarization profile of the source emission are 
often given in the local rest frame of the source. In Sec¬ 
tion 12 51 we use Eqs. (jH) to initialize the polarization 
propagator at the end of backward ray-tracing. 


2.3. Strong Lensing by Kerr Black Holes 

X-ray emission from an accretion disk is strongly lensed 
by the central Kerr black hole as it leaves the source to 
arrive at a distant observer. A lensed accretion disk ap¬ 
pears very different from a unlensed one. The frequency 
of the source emission is gravitationally and Doppler blue 
or redshifted depending on the position of the emitter, 
the 4-velocity of the source, and the distance and incli¬ 
nation angle of the observer. Consequently, the intensity 
prohle of a lensed accretion disk can be very different 
from an unlensed one, and the observed flux can also be 
different from the unlensed case. 

By definition, the monochromatic flux measured 
by an observer O very far from the black hole is (IMihalas I 

[T^ 

= j) = J I^^dLlo, (11) 

where dfto is the differential solid angle measured at the 
observer, and p is the angle cosine between the light ray 
and the normal of the detector window. We can safely 


where is the source frequency, g is the redshift factor of 
Eq. ([5|), {x°‘,Pa) is the photon’s position and momentum 
at the emitter found by backward ray-tracing (see Sec¬ 
tion 133) . To obtain Eq. (031) we have simply used the 
fact that Ivlv^ is conserved along each geodesic. The 
strong lensing magnification of the specific flux is then 

flensed 

P ^unlensed ' 

2.4. Backward Ray-tracing 

Since the emission region can be as close as a few gravi¬ 
tational radii from the black hole, the standard linearized 
gravitational lensing theory is not valid. Our strong lens¬ 
ing analysis in Kerr spacetime is based on backward ray¬ 
tracing. Given a source profile K(x°‘,p^), the flux inte¬ 
gral in Eq. m can be easily computed either analyti¬ 
cally or numerically. However, it is not a trivial job to 
evaluate the lensed flux integral Eq. (USD in which the 
integral is over the solid angle at the observer, but the 
integrand is evaluated at the source. We do this through 
backward ray-tracing. We choose a thin pencil beam 
(a sharp pentahedron) focused on the observer (the pig¬ 
ment core toward the black hole), and divide the solid 
angle space of this beam (large enough to contain all disk 
emission arriving at the observer) into a uniform grid of 
pixels. Through each pixel, we shoot one ray backward 
from the observer to the source. To compute the image 
area, we need only count the number of pixels whose light 
rays end at the accretion disk. To compute the observed 
flux, we weight pixels by the integrands in Eq. (1131) . In 
this way, our backward ray-tracing algorithm takes into 
account gravitational light deflection, Doppler or gravi¬ 
tational redshift, and area distortion simultaneously. 

2.5. Polarization and Gravitational Faraday Rotation 

Because the polarization at the source is known and 
needs to be computed at the observer, most current po- 
larization propagators are b ased on forward raytracing 
(|Schnittman fc Krolik l[^09t) . Our polarization integra¬ 
tor is based on backward ray-tracing. We assume that 
the polarization and limb darkening of the source emis¬ 
sion a re given as in the classical work of iChandrasekharl 
(|1960ll where the radiation transfer equation with po¬ 
larization was solved assuming a scattering-dominated 


























4 


semi-infinite atmosphere. Other initial polarization pro¬ 
files, such as optically-thin, inverse Compton scattering- 
dominated atmosphere, can also be implemented. 

At the end of the backward ray-tracing of each ray, 
we know the landing point on the accretion disk of a 
geodesic and its 4-momentum p“. First, we compute the 
photon’s 4 momentum components in the comoving 
frame. We then compute the angle cosine, /i, between the 
photon wave vector and the disk normal measured by the 
comoving observer. The angular dependent part of the 
intensity profile, w{^) (see Eq. (I2T|) in Section |4]), and 
the degree of polarization, Ssnur ri^lu), can be compute d 
by interpolating Table XXIV of iChandrasekhaFI (I1960D . 
We then compute the photon’s polarization vector E in 
the local rest frame using the fact that E is parallel to 
the disk plane and normal to the photon’s wave vector 
(£'-p = 0), i.e.. 


= 


VWV+W)‘- 




(15) 


Next we use Eq. ([S]) to compute the coordinate compo¬ 
nents of the polarization vector E. We then solve the 
parallel propagation equation of the polarization vector 
forward along the photon geodesic toward the distant 
observer using 


dE<^ 

~dr 


-TiyE^ 


(16) 


where F^^ are the Christoffel symbols. 

At the end of the forward polarization propagating 
process, we are back to the observer with rotated po¬ 
larization vector Eobs- Here we first compute the tetrad 
components of Eobs using Eq. ©■ We then work out the 
Stokes parameters corresponding to this ray. 


f & 

2y = arctan —^ 

V 

Q = cos 2x 
t7 = 5/y^ sin 2x (17) 

where x is the polarization angle, S = (^source is the de¬ 
gree of polarization which is conserved along a photon 
geodesic, and is the observed specific intensity at fre¬ 
quency r'o, which is related to the source intensity by 
= g^Ive ■ To compute the integrated degree of polar¬ 
ization and polarization angle, we add Q, U, and from 
each image pixel, obtaining Qtotai, Htotai, and Fi¬ 

nally we use Eq. (HZ} again to compute Xtotai and (5totai- 
Detailed algorithms are given in the next section. 


3. GUI AND NUMERICAL ALGORITHMS OF KERTAP 

3.1. Graphic User Interface 

The MATLAB version of KERTAP contains a GUI (see 
Figured]). It asks the user for the black hole parameters, 
the resolution desired, and the number of parallel work¬ 
ers for parallel computing. The GUI contains a status 
panel indicating the status of the code: waiting for in¬ 
put, running, or finished. If the code is run successfully, 
the GUI will generate the redshifted image of the accre¬ 
tion disk, and output important parameters such as the 
flux magnification and the degree and angle of the ob¬ 
served polarization. The performance of the code, such 


as the wall clock time used for ray-tracing, will also be 
recorded. 

3.2. Initialization Algorithm For Backward Ray-tracing 

Given an observer distance robs, and the dimensions 
of the image rectangle (orthogonal to the line of sight at 
the source redshift), Xsize and j/size, the angular size of the 
cone focusing at the observer can be inferred. Given the 
resolutions in the x and y directions, the pixel size can 
be computed as dxg^id = assize /^Tx and dygrid = 2/size/%- 
For a pixel labeled by {ix, iy) the photon’s starting 8D 
phase-space coordinates are computed in Algorithm 1. 


Algorithm 1 Initialization of one pixel (ix,iy) of the 
grid. are elements of the 3D photon momentum 

vector with respect to the ZAMO. Function p3_to_p4-() 
will normalize it assuming p* = 1, and then convert the 
ZAMO tetrad components to covariant coordinate com¬ 
ponents. We choose fobs = 0 and <j)ohs = 0 without loss 
of generality. 

Xq = (0, robs? ^obs» 0) ^ observer’s 4-coordinates 

p'^ = robs ^ must be positive 

= dfCgrid X ix 
p® = dpgrid X iy 

p4 = p3_£o_p4(p’",p^,p'^) o p4 is covariant 

Yq = (Xq,p4) > 8-D starting point for RT 


3.3. Ray-Tracing Using MATLAB ODE45 

For each ray arriving at the observer, we solve the 
geodesic equations using MATLAB ODE solver odefB. 
The syntax of ode45 is 

Ye,/e] = ode45(@F', grange, Tq, OpttOUS, ...), 

(18) 

where Y'{f^Y) is an 8-D column vector containing the 
right hand side (RHS) of Eq. dSj), Yq is the 8-D initial 
condition for ode45, and grange is the range of the in¬ 
dependent variable (affine parameter) where we want to 
solve the ODEs. A subtle point of backward ray-tracing 
is that we do not shoot rays toward the accretion disk 
(say, take p’' < 0), instead we must have p*" > 0 (toward 
the observer). This point is important because the Kerr 
metric is not symmetric with respect to the reflection of 
the time coordinate (time reversal will change the direc¬ 
tion of the black hole angular momentum). The back¬ 
ward ray-tracing is realized by properly choosing grange, 
i.e.,we define grange = [Co,6] where K(Co) = To, and 
y(Ci) corresponds to the light emitting event (Ci < Co, 
and < to). In Eq. (fT8l) options is a structure containing 
optional parameters that change the default integration 
properties, such as the relative and absolute error toler¬ 
ance, and events to stop the integration. For example, 
suppose we want to stop the ray-tracer when the geodesic 
is very close to the event horizon, say, r < (1 -|- 5)r_|_ for 
some small i5, we define a black hole hitting event as 

flagBH = (r < (1 + (5)r+) - 1. (19) 

Since we want to know when the geodesic hits the equa¬ 
torial plane, we define another event 

flagEquator = COS 9. (20) 

Then ode45 will record the point where flagen = 0 or 
flagEquator = 0 and stop the integration if the user re- 
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0^0 

Untitled 1 


1X1 KERTAP 


sRedshift Image of Accretion Disk= 



I 







= Performance!- 1 

Rate (rays/sec): | 105.105 | 

Wall Clock (hour); | 2,64814 1 


Finished Successfully 




Error_Q_max; ] 2.0136Sa-07 


: Paramete rs= 


spin a: 

inclination (deg): 
r_obs (unit M); 
r_disk (unit M); 
half X size: 
half y size: 
nx (half): 
ny (half); 
n umb_workers; 
photon index; 
radial index n; 
RUN 


— Output- 


g_min: 

g_max: 

Area Magnification: 
Flux Magnification; 
delta (polarization); 
chi (in degree); 


0,0321767 


0.0221874 


-6 56797 


Figure 1. GUI of KERTAP. Input parameters can be entered from the “Parameters” panel. KERTAP will run after the user pushes the 
“RUN” button. The status and performance of the program is shown in the “Status” and “Performance” panel, respectively. Important 
outputs are given in the “Output” panel. The redshifted image of a lensed accretion disk will be plotted after a successful run. Pushing 
the “Reset” button will clear the outputs and restore the default inputs. A parallel computing job can be submitted to a remote cluster 
directly from the GUI on the user’s desktop. 

quests to (as we did). These two events are defined in 
disk_events.m. In the left hand side of Eq. (US, the in¬ 
dex If, tells you which event happens, (^e,le) stores the 
point where the event happens, and {^,Y) stores the re¬ 
sults of the integration up to the event point (i.e.,the 
photon geodesic). Based on ode45, the ray-tracing part 
of KERTAP is easily realized in Algorithm 2. 


Algorithm 2 Ray-tracing Algorithm. MATLAB ode45 
is based on Runge-Kutta45. 

for ix = —nx : nx do 
for iy = —ny : ny do 

Initialize pixel (ix,iy) giving Yq See Algorithm 1 

[^, y, ...] = ode45(@y^(^), To, options, ...) 
if ray goes to infinity then 
go to next ray 

else if hit on the horizon then 
^hole “ ^hole 1 
go to next ray 

else if hit on the equatorial plane then 
if hit on the accretion disk then 
^disk ~ ^disk “1“ 1 

compute redshift, write (x“,p 5 ), ... 
go to next ray 
end if 
end if 
end for 
end for 


Algorithm 3 Polarization Propagating Algorithm. In 
the following, polarization () computes the polarization 
degree Ssource and polarization vector on the ac¬ 

cretion disk using {x°‘^p°‘) source obtained from backward 
ray-tracing. 5 source is interpolated from Table XXIV of 
iChandrasekhar I (jlQfiOf l. The polarization vector E°‘ is 
propagated forward from the source to the observer us¬ 
ing Eq. (HU). polar_obs() computes the polarization angle 
X at the observer, dobs = Ssource since the degree of po¬ 
larization is conserved along geodesics. 

for ix = —nx : nx do 
for iy = —ny : ny do 

[^,y, Ye] = ode45{@Y ',...) > Algorithm 2 

if a good shot then 

[Ssource, Esource] = polarization{Ye, M, a) 

Ecurr — Esource 
ripoint = length{^) 

for j — 0, f^point 2 do t> Esource ^ -^o6s 

for a = 1 : 4 do 

= - Etc=l l^curr) 

rpa _ Tpa \ j TPa 

Enext ~ Ecurr + dt, 
end for 
Picurr — Enext 
end for 
Eobs ~ Ecurr 

X = polar.obs{Eobs) 

end if 
end for 
end for 


3.4. Polarization Propagator 

At the end of the backward ray-tracing of each ray, 
the polarization propagator is initialized, and propagated 
forward to the observer, see Algorithm 3. 


3.5. Post Processing 

The results of ray-tracing are stored in different data 
structures. Eor example, Image_G and Image_Y store 
the redshifts {g = 0 for a missed shot) and 8-D phase 
space coordinates {x°‘,pb) of the target points on the ac- 
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cretion disk for all successful shots. Image_P stores the 
polarization data {w, S, x) where w is the limb dark ening 
factor from Table XXIV of I ChandrasekhaFI (|196ClD . The 
strong lensing amplifications of the observed flux and 
image area are computed in Algorithm 4. The observed 
degree and angle of polarization are computed using Al¬ 
gorithm 5. Thanks to the data visualization tool intrinsic 
to MATLAB and Python, the image of the lensed disk 
can be easily generated using functions such as imagesc() 
for MATLAB and imshow() for Python, and the gravi¬ 
tational Faraday rotated polarization vector field can be 
conveniently visualized using the function quiver() which 
is similarly defined for the two languages. The redshifted 
image of an accretion disk will be plotted in the GUI af¬ 
ter a successful run, and another two image windows (one 
for the same redshift image, the other for the intensity 
profile and polarization visualization) will pop up for the 
user to edit or save the figures in their desired formats 
(eps, jpg, etc). 


Algorithm 4 Strong lensing algorithm. dA and did are 
the differential image area and solid angle for one pixel. 
The unlensed image area and flux are easy to compute, 
and are assumed to be known. 

jplensed _ q. j^lensed _ q 

for ix = —nx : nx do 
for iy = —ny : ny do 
if a good shot then 

■'^lensed — -^lensed ^ 

t> Initialization 

g = ImagejG{ix^iy) 

> read the redshift g 

Y = ImageJY(ix,iy) 

> read 8-D 

dF = g^Pjg{Y) 

p'lensed _ jplensed _|_ 

end if 
end for 
end for 

j^lensed j^lensed y dA 

jplensed jplensed ^ 

> Us = Uojg 

,, _ Alenaed / Aunlensed 

parea — AA / Al 

> Area amplification 

_ jplensed jjpunlensed 

0 Flux amplification 


Algorithm 5 Polarization synthesizing algorithm, dtotai 
and xtotai are the integrated degree and angle of polar¬ 
ization. 

^total — 0) Qtotal — Oj Utotal — 0 

for ix = —nx : nx do 
for iy = —ny : ny do 
if IsImagePoint then 

0 Initialization 

g = Image_G{ix,iy) 

> read the redshift g 

r = ImageJY{ix, iy, 2) 

[w, S, x] = Image j^{ix, iy) 

> read r coordinate 

dl = wg^^^/r'^ 
dQ = S X dl X cos(2x) 
dU = S X dl X sin(2x) 

Qtotal ~ Qtotal ^Q 
^total — ^total “1“ 

Itotal — Itotal H“ (^I 
end if 
end for 
end for 

Xtotai — 0.5 tan {Utotal/Qtotal) 

^total — ^total/{Itotal ^'^^{‘^Xtotal)) 

> See Eq. ll23ll 


3.6. Python Version of KERTAP 

Despite the increased use of MATLAB in academia, for 
users who do not have access to a MATLAB cluster, we 
also provide a Python version of KERTAP. The Python 


code uses exactly the same algorithms as the MATLAB 
code except that it does not have a GUI. Furthermore, for 
the MATLAB code we used the ODE solver, odefS pro¬ 
vided by MATLAB (already highly optimized and with 
stringent error control), whereas for the Python code we 
created a 5*^ order Runge-Kutta solver using the Gash- 
Karp methods. The MATLAB ode45 provides a very 
convenient ‘event’ capturing mechanism to stop the ray¬ 
tracing at user defined events (such as when a ray crosses 
the event horizon, or punches through accretion disk, see 
Sec. ESI), whereas for the Python code, these events are 
coded into the ODE solver by hand. Consequently, to 
extend the code to include other events (e.g., higher or¬ 
der Kerr-lensing images) is slightly more difficult for the 
Python code than for MATLAB. 

The MATLAB version of KERTAP was parallelized us¬ 
ing the Parallel Computing Toolbox and the MATLAB 
Distributed Computing Server (MDGS). A very useful 
feature of the MDGS server is that it allows parallel jobs 
to be submitted from the user’s desktop directly to a dis¬ 
tant parallel computing clusters, i.e., the user provides 
the desired parameters through the GUI of KERTAP, 
and KERTAP automatically creates a job submission 
script and submits the job to the remote cluster deal¬ 
ing with popular job schedulers such as MOAB, SLURM, 
etc. The Python version of KERTAP was parallelized us¬ 
ing the free mpiJ^py package which supports a full imple¬ 
mentation of Message Passing Interface (MPI). The pure 
Python version of KERTAP is slower than the MATLAB 
version by a factor of ^2. We speed up the Python code 
by adding C-type static type declarations to the pure 
Python code using Cython which results in a ~100 times 
performance gain, see Fig. [T] We provide our users with 
a pure Python version of KERTAP, a Cython version, 
and a parallel Cython version. To run the Cython ver¬ 
sion of KERTAP, the Cython source file ‘mod.pyx’ which 
contains most of the routines needs to be compiled into a 
Python extension module ‘cymod.so’ before it can be im¬ 
ported as a Python package. To run the Cython code on 
a cluster in parallel, a MPI compiler (e.g., gnu-openmpi) 
is needed. Details of these implementations are included 
in the documentation to KERTAP. 


4. AN EXAMPLE APPLICATION 

We assume that the X-ray corona is a thin layer imme¬ 
diately above the AGN accretion disk, and moves with 
Keplerian flows (see Eqs. © and ©). We take rdisk = 
20 rg according to constraints obtained from quasar mi- 
crolensing observations (Kochanek 2004; Blackburne et 
al. 2006, 2014, 2015; Morgan et al. 2008, 2012; Chartas et 
al. 2009; Dai et al. 2010; Mosquera et al. 2013; MacLeod 
et al. 2015), and use nsco as the inner cutoff. The image 
window for backward ray-tracing is of size 50 x 50 rg. 
X-ray coronae of t his type are motivated by t he so- called 
sandwich model (jHaardt fc Maras^ll99lL 1199311. and 
are o ften used in the literature ( Ruszkowski fc Fabian I 
[200(1 . For a corona geometry such as a sphere, or a 
disk above the accretion disk, see Chen et al. (2013a). 
Because quasar X-ray emission is observed to follow a 
power law, we assume the following form for the source 
intensity profile. 


/^(z/,^, r) oc 


1 w(m) 


( 21 ) 
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where fj, is the cosine between the photon wave-vector 
and the upward disk normal measured by the comov¬ 
ing observer, and w{^) is the angular-dependence of the 
intensity profile (we have assumed azimuthal symmetry 
for simplicity). To be more specific, we take from 
[Chandrasekhar I l)196fl[ l (classical results from solving the 
radiation transfer equation with polarization assuming a 
scattering-dominated semi-infinite atmosphere). The pa¬ 
rameter r is called the photon index, and n specifies the 
radial steepness of the profile. We take T = 2.0 (e.g., 
Chen et al. 2012), and n = 3 (radially steep) or 0 (ra¬ 
dially flat). We test for two black hole spins, a = 0.998 
(Thorne 1974) or 0 (Schwarzschild black hole). We take 
the observer distance as robs = and test for differ¬ 

ent inclination angles (from nearly face on to nearly edge 
on). 

For the assumed simple geometry and intensity profile, 
the unlensed specihc flux at frequency Vo can be analyt¬ 
ically integrated using Eq. (HID 


F, 


unlensed 


27r/iobs^/^(A^obs) 1 

,.r-i 


' obs 


dr 

^n —1 


2TTflohsWifiohs) 1 

(«-2)rL 


nsco 

1 

disco 


j.n-2 

' disk -I 


( 22 ) 


where /iobs = cos(0obs)- We have dropped the unimpor¬ 
tant constant in the definition of 1^. For n = 2, the {n—2) 
factor on the denominator is absent, and the expression 
within the square brackets is replaced by logrbisk/usco- 
From Eq. (d the lensed flux is 


F, 


lensed 




r-1 


dQ.n 




,T-1 

' obs 


C^^grid, (23) 


where the summation is over the accretion disk image, 
and dAgrid = dxg^id x dyg^id is the area of one image 
pixel. Note that 'w{fi) varies from pixel to pixel because 
of the general relativistic light bending, and the aberra¬ 
tion effect. The magnification p, = ^^ensed^^miiensed 

be computed using Eqs. (US) and (1^ immediately (note 
that the term ^cancels out). 

In Figure [2] we plot the lensed intensity profile and dis¬ 
tribution of polarization for the X-ray disk observed at 
inclination angle 75°. Figure |3] is similar to Fig. [21 but 
is for black spin a = 0.5 and we have included up to 4*^ 
order Kerr-lensing images. Kerr strong lensing changes 
both the shape and the intensity profile of the disk signif¬ 
icantly. The intensity is strongly concentrated in a small 
region (in dark red) near the black hole and to the left 
of the accretion disk (where the rotating source is ap¬ 
proaching the observer). As for polarization, the degree 
of polarization to the left of the black hole is on aver¬ 
age smaller than the classical value (~4.6% for 9 = 75°; 
Chandrasekhar 1960) because the angle between the pho¬ 
ton wave vector and the disk normal measured by the 
comoving observer is smaller than the inclination angle. 
Therefore, ^(/i) is smaller, see the black lines in the sec¬ 
ond row of Figure |T| for the inclination dependence of 
5. The contrary holds for the right hand side of the 
accretion disk. On the other hand, the direction of po¬ 


larization can change significantly because of the gravi¬ 
tational Faraday rotation. This effect is more significant 
for source regions closer to the black hole where the ob¬ 
served polarization angle x can be positive or negative 
(depending on the act ual source location) inst ead of hor¬ 
izontally aligned as in [Chandrasekhar I ([196r)tl . 

The strong lensing magnification of the image area and 
specific flux, and the integrated (or averaged) polariza¬ 
tion degree and angle are plotted in Figure HI We did 
not consider higher order images here since the accre¬ 
tion disk is believed by many to be optically thick. We 
plot the inclination dependence of //area, //flux, 5, and y 
for spin a = 0.998 and 0, and for steep and flat radial 
prohle (n = 3 or 0). The magnification increases with 
inclination angles and is more significant for steeper pro¬ 
files and for larger spins (smaller risco) both of which 
amount to more concentrated emission closer to the black 
hole. The fact that for small inclination angles (nearly 
face on), the observed degree of polarization is greater 
than the classical result is mainly caused by gravitational 
light bending ((1 = 0 for face on case, see Chandrasekhar 
1960). For moderate to high inclination angles, the emis¬ 
sion is strongly focused in a small region with a low de¬ 
gree of polarization, see Figure |2j Consequently, special 
and general relativity effects tend to reduce the observed 
degree of polarization confirming earlier results (Connors 
et al. 1980). 

Another interesting point is that for a fixed emission 
size (e.g., rdisk = 20 as used in this paper) with the 
intensity profile Eq. (12111 . the effects on polarization, 
i.e.jthe polarization profile, the integrated degree and 
the angle of polarization, are achromatic. This is because 
both the input s of the polarization pr opagator, i.e.,(5(/t) 
and w{y,) from [Chandrasekhar I (jl960D . and gravitational 
Faraday rotation are wavelength independent, and the 
factor cancels out when averaging over the im¬ 
age of the disk when computing Xtotai and iJtotai, see 
Eqs. (fTTI). (1221). and (1221). This apparent discrepancy 
with [Schnittman &: Krolikl ([200^ is caused by the fact 
that [Schnittman fc KroliO ^ ()2009[ ) have assumed thermal 
emission for X-rays emitted by the accretion disk around 
a stellar mass black hole, and therefore, the X-ray emis¬ 
sion of higher frequencies are emitted by regions closer 
to the black hole (the effective temperature of the disk is 
higher closer to the black hole), whereas we have assumed 
the same power-law oc at each radius for X-ray 

emission by AGN. If we assume that hard X-rays are 
emitted by a region smaller than that of th e soft X-rays 
and cl oser to the black hole as suggested bv [Chen et al. I 
unni), say, take rdisk = 10 then, even under the same 
power-law assumption, the integrated degree and angle 
of polarization should differ from current results, i.e., the 
results will become chromatic. 

The accuracy of KERTAP is very good. We have used 
Carter’s constant as an independent check of accuracy. 
At the end of the long ray-tracing (robs = 10®rg), the 
relative error of Carter’s constant is never worse than 
~ 10“^ for all cases tried by the authors. We also checked 
the accuracy of KERTAP using the norm of the photon’s 
4-momentum vector and the polarization vector £1“, see 
Fig. [5] for results of a typical ray. The results of the 
MATLAB code are highly consistent with the Python 
code, e.g., the norm of the difference between the red- 



































Figure 2. Intensity and polarization profile of a X-ray disk of ra¬ 
dius r = 20 strongly lensed by a Kerr black hole of spin a = 0.998 
observed at 0 = 75°. The colorbar is given in log scale. Without 
the strong lensing effect, the image is an ellipse (with a hole), and 
the intensity profile is symmetric with respect to the y-axis. The 
lensed disk is warped upward, and the intensity is concentrated in 
a small region on the left hand side of the disk (approaching the 
observer). The polarization should be horizontally aligned and con¬ 
stant in value ignoring spacetime curvature (5 = 4.6% for 0 = 75°, 
see Chandrasekhar 1960). The strong lensing of the Kerr black 
hole reduces the degree of polarization significantly (-^1.6%). The 
orientation of the polarization can also be rotated near the black 
hole. 

shift images generated using the Python and MATLAB 
code is of order and there were no rays which hit 

on the accretion disk or enter the horizon in one case 
that did not hit on the disk or entering the horizon in 
the other. A convergence test is shown in Table [T] and 
Figure [6] The convergence of the KERTAP is pretty 
good, e.g., the amplification to image area and specihc 
flux, and the polarization degree and angle all converge 
at a grid of moderate size 500 x 500. In Fig. Owe test 
the scalability of KERTAP. Both the MATLAB (par¬ 
allelized using MODS) and Python (parallelized using 
mpi4py) code scale roughly linearly with the number of 
CPUs. Consequently, the performance of KERTAP can 
be significantly improved by parallel computing. How¬ 
ever, the performance improvement of the MATLAB 
code through parallel computing was restricted by the 
number of MDCS licenses available in a user’s MATLAB 
cluster, whereas the Python code is only restricted by the 
size of the cluster (number of CPUs available for the MPI 
job). It is important to note that since MATLAB20I4b, 
the Parallel Computing Toolbox allows up to 512 parallel 
workers (this number is not constrained by the number 
of MDCS licenses available). Consequently, significant 
speedup can still be achieved by users without a MAT¬ 
LAB cluster but with a single node with many cores. 

5. DISCUSSION 

We have developed KERTAP, the first ray-tracing code 
using interpreted languages MATLAB and Python, for 
studying strong gravitational lensing in Kerr spacetime 
including polarization^ The key ingredients of KER¬ 
TAP are: a GUI (for the MATLAB version), backward 
ray-tracing realized by 5*^ order Runge-Kutta method, a 
polarization propagator and synthesizer, and the formal¬ 
ism and algorithms for computing strong lensing effects 
of a Kerr black hole centered on an accretion disk. The 
backward ray-tracing part of KERTAP most resembles 



Figure 3. Intensity and polarization profile of a X-ray disk of 
radius r = 10 strongly lensed by a Kerr black hole of spin a = 0.5 
observed at 0 = 75°. We include up to 4^^ order Kerr-lensing 
images. 

that of Schnitt man &: Bertsch i nger ( 2004) and Vincent 
et al. (2011). iVincent et al. I (|2011h focuses on radia¬ 
tion transfer while this paper focuses on observable ef¬ 
fects of Kerr strong lensing and polarization. The in¬ 
dependence of the geodesic equations from the polar¬ 
ization propagation equations makes it possible to solve 
the geodesic equations (backward) first, then propagate 
the polarization vector forward. The idea of backward 
ray-tracing followed by forward polarization propagating 
is new, to the best of our knowledge. Polarization and 
gravitational Faraday rotation have been studied exten¬ 
sively in the literature (e.g., Agol 1997; Agol & Kro- 
lik 2000; Schnittman & Krolik 2009). However, pub¬ 
lic code explicitly dealing with polarization is not yet 
available. Published papers investigating the X-ray po¬ 
larization of accretion disks are basically all doing for¬ 
ward ray-tracing ( with or without Monte-Carlo r adiation 
transfer; see, e.g., ISchnittman fc: Krolik I (|2009ll l. Since 
X-rays are emitted by sources very close to the central 
black hole and moving with highly relativistic flows, in 
order to compute the observable property measured by 
a distant observer at a given inclination angle, at each 
point on the accretion disk and for each direction, many 
forward rays have to be traced (because no one knows 
which ray will go where). This increases the forward 
ray-tracing time significantly, and the results are usu¬ 
ally noisier than tho se from backward ray-tracing (e.g., 
compare Figure 1 of ISchnittman fc Krolik I (1200911 with 
Figure [2] of this paper). 

The idea of performing very computationally intensive 
jobs such as ray-tracing in Kerr spacetime using inter¬ 
preted languages such as MATLAB or Python might 
sound odd at the first. It is hard to expect a dynamic 
language such as Python to compete in performance with 
FORTRAN or C-I--1-3 In fact, the single core (MATLAB 
or pure Python) version of KERTAP performs poorly 
compared to code written in FORTRAN and C-|—1-, and 
to GPU-accelerated CUBA code (see, e.g., Dexter & Agol 


^ The code is available 

https: //bitbucket.org/binchenl4/kertap 


at ^ This does not imply that the back and forth ray-tracing algo¬ 

rithms presented in this paper are intrinsically slow. 
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Figure 4. Inclination dependence of Kerr strong lensing and polarization. We test for two black hole spins a = 0.998 and 0 (the first and 
second column), and two radial profiles n = 3 (steep) and n = 0 (fiat). The first row plots the strong lensing magnification of the image 
area and total flux. The secon d and third rows plot o bserved degree and angle of polarization, respectively. The black lines in the second 
and the third row are classical IChandrasekharl il960l) results. 



Figure 5. Accuracy of KERTAP. A typical ray is shot backward 
toward the black hole, hits on the accretion disk, after which the 
polarization vector is forward propagated to the observer. The 
dot-dashed blue curve shows the relative error in Carter’s con¬ 
stant. The dashed magenta curve shows the norm of the photon 
4-momentum vector. The solid red line shows error in the norm 
of the forward parallel-propagated (space-like) photon polarization 
vector. 


2009; Chan et al. 2013). However, this does not mean 
that interpreted languages such as Python or MATLAB 
will/should not be used in CPU intensive astronomi¬ 
cal computations such as Kerr ray-tracing. The perfor¬ 
mance of these languages has been improved significantly 
over the past decade, in particular, through the devel¬ 
opment of application programming interfaces (APIs) 
with low-level languages such as C/FORTRAN (e.g., the 
CPython API), multi-CPU parallel computing (mpidpy 
for Python; MDCS for MATLAB), and by supporting 
GPUs. For example, both Python and MATLAB sup¬ 
port GPU programming, e.g., PyCUDi^ for Python. 
The performance of the pure Python version of KER¬ 
TAP was improved by two orders of magnitude using 
Cython. The disadvantage in speed is compensated by 
parallel computing. In MATLAB, this is achieved by 
using the Parallel Computing Toolbox on local machines 
with multiple cores, or by using the Distributed Comput¬ 
ing Server over a cluster. For Python, the Cython version 
of KERTAP supports full implementation of MPI. For 
example, we have shown in this paper that it is feasible 
to study strong lensing of X-ray emission of accretion 
disks around Kerr black hole using MATLAB/Python. 
The strong lensing analysis of a X-ray emitting disk of 

® http://mathema.tician.de/software/pycuda/ 
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Figure 6. Convergence test of KERTAP. We plot the resolu¬ 
tion dependence of the maximum and minimum redshift, ^max and 
5'mirn the strong lensing magnification of the image area and spe¬ 
cific fiux, /i-area and pflux? and the observed degree and angle of 
polarization. The data are tabulated in Table U For an X-ray 
emitting disk of size 20rp, a 500 x 500 grid for an image rectan¬ 
gle of size bOvg X bOvg is accurate enough for the Strong lensing 
calculations presented in this paper. A grid of this size takes ^^45 
minutes to run on a single computer with 8 MATLAB workers. 
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Figure 7. Scalability of the parallel Python and MATLAB code. 
The test was run for the example in Fig. [2] with assize = 25 
Vsize = 15 Tg, and resolution 800 x 480. The rays were shot from r = 
10® Vg from the black hole, with a stringent error control (relative 
error e = 10“^^ for the RK45 ode solver). The Python code was 
parallelized using the ‘mpi4py’ package whereas the MATLAB code 
was parallelized using the MDCS. Both codes show good scalability. 
The Python code was sped up by two orders of magnitude using 
Cython. The curve is steeper for the MATLAB code when the 
number of CPUs is small, this is simply caused by the fact that 
among the N-workers allocated for the parallel job, one is selected 
as the master worker, and only N-1 workers are doing the ray¬ 
tracing. The number of parallel MATLAB workers is restricted 
by the number of MDCS licenses available on a MATLAB cluster, 
whereas the number of Python workers is only restricted by the 
size of the cluster. 


radius 20 can be done with high accuracy in ~40 min- 
utes by running the Matlab version of KERTAP on a 
desktop computer with multiple cores (< 1 minute using 
the Cython code). The computing time can be further 
reduced on a high performance computing cluster given 
the massively parallel nature of ray-tracing and the good 
scalability of KERTAP (see Pig.0. The MATLAB GUI 
of KERTAP allows a ray-tracing job to be created on a 
local machine (e.g., the user’s laptop), submitted to a 
remote high-performance computing cluster with the re¬ 
sults sent back to the user upon a successful run, the most 
desired form of interactive high performance computing^ 
If the reader desires, the strong lensing and polarization 
algorithms can be easily realized in FORTRAN or C-I--I-. 

The advantage of coding using high-level languages 
such as Python and MATLAB is multifold. First, the 
object-oriented nature of MATLAB makes it very easy 
to build graphic user interfaces. For example, the GUI 
as shown in Figure [T] is written using the MATLAB 
GUI building tool “guide.” Secondly, it is much easier to 
write code using MATLAB or Python than FORTRAN 
or G/G-I--I-. Programming in MATLAB/Python greatly 
shortens the whole cycle of a scientific project, i.e., from 
prototyping and debugging the code, to production runs, 
to visualization of the scientific results and generating 
high quality plots. For example, this whole paper can 
be done using a single language (either MATLAB or 
Python). Our code is very short. The kernel of the 
MATLAB version of KERTAP is only about 600 lines. 
The kernel of the Python version of KERTAP is about 
200 lines longer (including the Runge-Kutta solver). The 
reader should experience no difficulty in understanding 
our code after reading this paper. Thirdly, since MAT¬ 
LAB was designed to be a data analysis and visualization 
tool (and thanks to free Python packages such as ‘mat- 
plotlib’), the post-processing of ray-tracing data using 
MATLAB or Python is almost trivial. For example, the 
piece of code generating the redshifted image of the ac¬ 
cretion disk (see Figure [T]) is only a few lines. Users do 
not have to write separate codes visualizing their results 
using software such as IDL or Mathematica. Our code 
would be very useful for scientists new to the Kerr geom¬ 
etry or in need of polarization calculations, or to those 
people who need a second code to compare with their 
own. KERTAP is the hrst public code explicitly dealing 
with strong lensing of X-ray polarization. For less com¬ 
plicated computing jobs, our MATLAB/Python code can 
be used directly or after some modification. For compli¬ 
cated computing jobs, KERTAP might still be useful if 
the user’s institute has a parallel computing cluster since 
KERTAP is fully parallelized. We conclude that it is 
possible to perform complex numerical-relativity-related 
computations using interpreted languages such as MAT¬ 
LAB and Python. 
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Table 1 

Convergence and performance of KERTAP. We run KERTAP on a single 
computer with eight MATLAB parallel workers. Accretion disk size 
^disk = 20 Tg. The size of the image window is 50 x 50 r^. Black hole spin 
a = 0.998. Photon index P = 2.0 and radial profile n = 3. Inclination angle 
0 = 75°. Observer distance robs = 10^ 'f'g- 


Resolution 

P'area/P'flux 

5min/Pmax 

& 

X 

(deg) 

Wall time 
(hr) 

Speed 
(ray s”!) 

20 X 20 

1.563/1.742 

0.257/1.280 

0.0230 

-8.473 

0.004 

27.33 

50 X 50 

1.568/1.587 

0.051/1.334 

0.0167 

-6.728 

0.011 

66.15 

100 X 100 

1.567/1.712 

0.039/1.356 

0.0159 

-8.017 

0.034 

83.65 

200 X 200 

1.563/1.723 

0.037/1.357 

0.0158 

-8.052 

0.120 

93.07 

500 X 500 

1.563/1.726 

0.038/1.360 

0.0158 

-8.005 

0.713 

97.75 

1000 X 1000 

1.563/1.726 

0.036/1.360 

0.0158 

-7.989 

2.881 

96.62 






